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SECURITY 


SUMMARY 


The  purpose  of  the  Trajectory  Reconstruction  and  Analysis  Method  (TRAM)  project 
is  to  develop  software  which  provides  the  Western  Space  and  Missile  Center 
(WSMC)  with  an  improved  Best  Estimate  of  Trajectory  (BET)  analysis  capability. 
Data  from  on-board  and  land  based  sensors  have  noise  content  of  varying  degrees. 
In  order  to  derive  information  from  this  noisy  data  optimally,  the  TRAM  program 
must  be  computational ly  quite  complex  and  require  significant  computer  time 
when  run  on  IBM  360/65  or  an  equivalent  machine.  Because  of  this,  algorithms 
which  economize  the  computation  have  been  developed. 

The  vehicle  trajectory  and  transition  matrix  computation  consume  an  appreciable 
portion  of  the  filtering  and  smoothing  process  of  TRAM  and,  consequently, 
efficient  methods  to  compute  these  quantities  were  needed  and  developed. 

However,  prior  to  implementing  these  efficient  algorithms  in  the  software  code 
a  'standard'  set  of  algorithms  that  do  not  degrade  the  accuracy  of  the  data 
and  preserve  the  accuracy  of  the  estimation  process  must  be  constructed.  The 
proposed  set  of  efficient  algorithms  are  then  compared  to  the  standards  to 
insure  that  the  efficient  methods  are  sufficiently  accurate. 

This  study  undertakes  to  establish  standards  in  two  areas.  The  first  is  the 
computation  of  the  powered  flight  trajectory  used  to  compute  the  nominal  tra¬ 
jectory,  the  nominal  radar  off-set  data  and  the  transition  matrices.  The 
second  area  in  which  standards  are  established  is  in  the  computation  of  the 
powered  flight  transition  matrices. 

Two  integrators  were  used  to  integrate  the  powered  flight  trajectory,  and  sub¬ 
sequently,  values  within  the  interval  of  integration  were  obtained  by  interpola¬ 
tion.  These  methods  took  into  account  a  trajectory  having  an  appreciable 
acceleration  component  due  to  the  thrust  of  powered  flight  and  a  trajectory 
reconstructed  from  quantized  acceleration  data  that  limits  its  precision. 


The  powered  flight  transition  matrices  are  subdivided  into  several  submatrices 
in  [2],  The  present  study  deals  only  with  the  matrices  designated  by  <j>  in 
Reference  [2]. 

The  character i stic  matrices,  which  are  used  in  the  computation  of  the  transition 
matrices,  tested  both  the  Runge-Kutta  and  Adams-Moul ton  integrators  for  integra¬ 
tion  accuracy.  The  Runge-Kutta  method  proved  markedly  better  and  met  the  accuracy 
criterion.  The  optimal  step  size  in  this  study  turned  out  to  be  as  large  as 
the  one  used  in  the  free  fall  portion  of  TRAM.  This  result  is  demonstrated 
mathematically  in  this  technical  report. 

By  integrating  the  equations  at  very  small  step  sizes,  one  may  obtain  an  estimate 
of  the  accuracy  of  the  foregoing  algorithms.  The  results  of  the  integrations 
at  larger  step  sizes  and  the  results  from  interpolation  were  compared  to  those 
resulting  from  integration  at  small  step  sizes.  This  comparison  thus  determines 
an  estimate  of  accuracy  for  both  integration  and  interpolation. 


Reference  [2],  B  \As,  R.  A,,  Trajectory  Reconstruction  and  Analysis  Methodology, 
Vandenberg  Air  Force  Base,  Performance  Analysis  Department,  Federal  Electric 
Corporation,  WTR  Division,  1978. 
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INTRODUCTION 


i.O 

This  report  is  an  extension  of  an  earlier  work,  Reference  [1],  which  was  con¬ 
cerned  with  the  efficient  and  accurate  calculation  of  the  state  and  transition 
matrices  for  the  TRAM  project,  Reference  [2].  The  earlier  work  was  concerned 
with  the  vacuum  freefall  portion  of  the  trajectory;  whereas,  this  work  is  con¬ 
cerned  with  the  powered  flight  portion. 

Because  of  the  demands  placed  on  the  computer  in  the  environment  in  which  a 
program  like  TRAM  must  be  executed,  there  is  always  a  great  need  to  make  the 
execution  as  efficient  as  possible.  In  addition,  the  numerical  algorithms  must 
be  sufficiently  precise  so  that  the  accuracy  of  the  instrumentation  and  the 
inherent  capability  of  the  estimation  method  are  not  degraded.  This  study  was 
undertaken  to  guarantee  that  result. 

Two  important  differences  between  powered  flight  and  freefall  are  that  the 
dynamics  of  the  vehicle  are  greater  in  powered  flight  because  of  the  thrust 
forces  and  that  information  about  its  acceleration  due  to  thrust  is  supplied 
by  PIGA  counts.  However,  the  PIGA  counts  are  quantized  so  that  their  accuracy 
is  limited  and  they  are  available  only  at  discrete  times.  This  is  in  contrast 
to  freefall  where  the  information  about  the  trajectory  supplied  by  the  differ¬ 
ential  equations  of  motion  is  always  continuous. 

Because  the  dynamics  are  greater  and  the  differential  equations  are  discontinu¬ 
ous,  it  is  clear  that  it  will  be  more  difficult  and  costly  to  obtain  the 
required  accuracy  in  powered  flight.  Therefore,  the  numerical  methods  developed 
for  powered  flight  must  be  able  to  handle  the  greater  dynamics  and  to  make  the 
optimum  use  of  the  data  available. 

The  transition  matrix  composition  differs  in  powered  flight  from  free  fall. 
Reference  [2].  It  is  necessary  to  include  terms  related  to  the  inertial 


Reference  [1],  Thompson,  G.  T. ,  Computation  of  State  Vector  and  Transition 
Matrix  for  TRAM,  Vandenberg  Air  Force  Base,  Federal  Electric  Corporation,  WTR 
Division,  1977. 

Reference  [2],  Brooks,  R.  A.,  Trajectory  Reconstruction  and  Analysis  Methodology, 
Vandenberg  Air  Force  Base,  Performance  Analysis  Department,  Federal  Electric 
Corporation,  WTR  Division,  1978. 
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measuring  unit  in  the  state  vector.  These  terms  are  static;  nevertheless, 
they  ao  affect  position  and  velocity  and  the  transition  matrix  must  reflect 
this.  Tin's  report  is  concerned  only  with  the  part  of  the  transition  matrix 
which  is  associated  with  the  dynamic  states  of  position  and  velocity.  The 
part  associated  with  the  IMU  static  states  will  be  treated  elsewhere. 
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2.0 


PROBLEM  STATEMENT 


2.1  The  Trajectory  Equations 

In  both  the  freefall  and  powered  flight  segments  of  the  trajectory  the  estima¬ 
tion  is  mechanized  around  a  nominal  trajectory  which  can  be  obtained  by  inte¬ 
grating  the  equations  of  motion.  The  nominal  trajectory  is  used  to  compute  a 
nominal  variation  offset,  the  transition  matrices,  and  it  is  used  to  compute 
the  total  state  vector. 

The  equations  of  motion  are  expressed  formally  as: 

P  =  V 

V  =  Ag(P)  +  Ar(P,  V)  +  Ar( t)  (2. 1) 

where  P  and  V  are  three  component  vectors  of  position  and  velocity  respective¬ 
ly;  A^  is  the  acceleration  due  to  gravity;  AR  is  the  apparent  acceleration  due 
to  the  rotation  of  the  earth,  and  A^  is  the  acceleration  due  to  thrust.  Gravity, 
as  is  indicated,  is  a  function  of  position.  The  acceleration  due  to  the  rota¬ 
tion  of  the  earth  is  a  function  of  both  position  and  velocity.  Thrust  acceler¬ 
ation  is  a  function  of  time  alone.  During  freefall  Ay  is  identically  zero. 

During  powered  flight  it  is  determined  from  the  PIGA  data  telemetered  every 
0.03  seconds.  The  accuracy  of  the  solution  of  equation  set  (2.1)  is  investi¬ 
gated  in  this  report. 

2.2  The  Transition  Matrices 

Let  x  be  a  six  component  vector  whose  first  three  components  are  position  and 
whose  last  three  components  are  velocity.  Then  equations  (2.1)  can  be  written 

x  =  f ( x ,  t) 

Further,  let 
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F(t) 


A  3  f 

3x^(x,  t) 

Then  the  transition  mat""'  satisfies  the  differential  equation 

<t>(t,  tQ)  =  F(t)«Kt,  tQ)  (2.2) 

where  <J»(t  ,  t  )  -  I 

and  I  is  the  identity  matrix. 

The  transition  matrix  will  have  an  order  which  is  equal  in  general  to  the  size 
of  the  state  vector.  However,  in  this  study  it  will  only  have  an  order  of  six, 
because  this  study  is  only  concerned  with  the  part  of  the  transition  matrix 
associated  with  the  dynamic  states  of  position  and  velocity.  The  accuracy 
with  which  equation  set  (2.2)  can  be  solved  is  investigated  in  this  report. 

2.3  Integration  Accuracy 

The  solution  of  the  differential  equation 

x  =  f(x,  t) 

x(0)  -  xq 

can  be  expressed 

*r  (t  j  -  xT(t)  +  r.  R(t)  +  f-T(  t ) 

where 

xc(tf  ,  the  computed  value  of  x, 

Xj(t)  is  the  true  value. 
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f.p(t)  is  the  round-off  error, 
and  r-|.(t)  is  the  truncation  error. 

<y(t)  results  from  the  computer's  use  of  an  approximation  to  the  derivative, 
results  because  of  the  finite  word  length  of  the  computer. 

Let  t:  -  r.p  +  t.y ,  then  the  following  graph  is  representati  ve  of  the  total  error 
(.  at  a  fixed  time  as  a  function  of  the  step  size  h. 


•r 


In  Figure  1  h  is  the  integration  step  size.  This  figure  indicates  that,  in 
general,  there  exists  an  optimum  h  at  which  the  total  error  has  a  minimum.  If 
the  step  size  is  increased  from  the  optimum,  then  the  total  error  increases 
because  Cy  increases,  even  though  decreases.  If  the  step  size  is  decreased 
from  the  optimum,  then  the  round-off  error  causes  the  total  error  to  increase 
even  though  the  truncation  error  decreases.  The  determination  of  an  optimum  h 
is  the  purpose  of  this  study.  Clearly,  it  is  desirable  to  use  as  large  a  value 
of  h  as  possible  for  which  the  total  error  is  within  allowable  bounds,  because 
larger  the  values  of  h  result  in  smaller  execution  times. 

2.4  Polynominal  Interpolation 

Both  the  trajectory  and  transition  matrices  will  be  needed  at  arbitrary  points. 
Normally,  and  specifically  for  the  cases  under  consideration  here,  the  integra¬ 
tion  of  the  differential  equations  is  much  more  expensive  than  evaluating  in¬ 
terpolating  polynomi nal s .  It  is  therefore  more  efficient  to  integrate  the 
respective  differential  equations  over  as  large  an  interval  as  possible,  and 
then  interpolate  for  values  within  these  intervals. 

Since,  after  integrating  the  differential  equations,  the  function  and  its 
derivative  are  readily  available  at  both  ends  of  the  interval,  spline  poly¬ 
nomials  are  obvious  choices  as  interpolators  because  the  spline  function  and 
its  derivatives  must  be  continuous  at  the  interval  mesh  points.  However, 
whatever  choice  is  made,  the  interpolator  must  contribute  errors  which  are 
significantly  less  then  the  error  of  the  integrator.  The  determination  of  the 
error  of  interpolation  is  one  purpose  of  this  study. 
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3.0 


ANALYSIS 


3.1  Analysis  Of  The  Equations  Of  Motion 

3.1.1  Formulation  Of  The  Equations  Of  Motion 

There  were  two  approaches  made  in  formulating  the  equations  of  motion.  The 
first  is  designated  as  Set  I.  It  is  that  set  given  in  Se;tion  2.0 

P  =  V 

V  =  Ag(P,  V)  +  Ar(P,  V)  +  Ap( t)  .  (3.0) 

The  Ag(P,  V)  term  is  a  simple  gravity  model  including  a  J  and  a  D  term.  A„(P,V) 
is  the  acceleration  present  because  the  equations  are  expressed  in  a  rotating 
coordinate  system.  It  includes  the  centrifugal  and  centripital  terms. 

Ap(t)  is  the  acceleration  due  to  thrust.  It  is  computed  from  the  PIGA  counts, 
counts  from  a  pendulous  integrating  gyroscopic  accelerometer.  Since  the  PIGA 
counts  give  the  change  of  velocity  over  the  minor  cycle,  all  that  needs  to  be 
done  to  determine  the  average  acceleration  over  the  minor  cycle  is  to  divide 
the  change  in  velocity  by  the  length  of  the  minor  cycle,  0.03  seconds.  The 
PIGA  counts  are  converted  from  counts  to  feet  per  second  by  multiplying  by  the 
quantization  factor  of  0.12,  and  the  results  are  transformed  from  the  instrumen¬ 
tation  coordinate  system  to  the  earth  centered  rotating  coordinate  system  in 
which  the  equations  of  motion  are  integrated. 

The  Ap  term  is  discontinuous,  and  therefore  adds  to  the  difficulty  of  integrating 
(3.0)  with  as  much  accuracy  as  desired,  so  an  alternate  approach  was  formulated. 
It  is  designated  as  Set  II. 

First  the  position  and  velocity  resulting  from  thrust  is  calculated. 


t 

Vp(t)  =  fQ  Ap(a)do 


(3.1) 
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t 

PF(t)  =  /0  VF(l)dT 

Then  primed  variables  are  introduced 

V1  J  V  -  Vp  (3.2) 

P'  t  P  ’  PF  • 

With  these  variables  the  differential  equations  of  motion  are  written  as: 

P'  =  V* 

V'  =  Ag(P'  ♦  PF)  ♦  Ar(P*  +  PF,  V'  ♦  VF)  .  (3.3) 

Finally,  the  total  position  and  velocity  are  computed  by 

V  =  V1  +  Vp  (3.4) 

P  =  P'  +  Pv  . 

From  (3.3)  it  is  seen  that  the  discontinuous  term  is  not  present  in  the  dif¬ 
ferential  equations.  This  formulation  will  aid  in  the  determination  of  the 
solution  accurately. 

The  terms  Vp  and  Pp  are  calculated  directly  from  the  PIGA  data,  and  so  the 
method  fo’-  doing  this  needs  to  be  set  forth  explicitly. 

Since  the  PIGA  counts  represent  changes  in  velocity,  it  is  natural  to  represent 
the  counts  for  a  i  nor  cycle  as  AV..  Therefore, 

W  ;  .\  AVi 
1  =  1 

because  Vp(0)  =  0. 
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The  position  due  to  thrust  is  approximated  by  integrating  velocity  by  the 
trapezoidal  rule. 

PF(tn)  =  PF(tn.1)  +  MVF(tn)  ♦  vF(tn_1))( . 03) 

Pp(0)  -  o  . 

The  integration  step  size  is  the  minor  cycle  time  of  0.03  seconds. 

Ap ( t n )  is  computed  by 
Ap(tn)  -  AVn/0.03  . 

By  the  use  of  the  above  formulas  Ap(tn),  Vp(tn)  and  Pp ( t p )  n  -  0,  1,  2  .... 
can  be  computed.  The  values  at  any  time  can  be  found  by  interpolating,  as 
fol lows: 

Ap(t)  =  AF(t.)  ,  t.  <  t  <  ti+1 

VF(t)  =  VF(t  - )  ♦  (t  -  t.)AF(t.)  t.  ^  t  -  t.  +  1 

PF(t)  -  PF(t.)  Mt  -  ti)vF(t.)  +  s»(t  -  ti)2AF(ti)  t.  i  t  ^  ti+1 

The  results  of  the  interpolation  are  multiplied  by  the  scale  factor  to  change 
them  from  counts  to  feet  per  seconds,  and  then  they  are  transformed  from  the 
instrumentation  coordinate  system  to  the  earth  centered  rotating  system. 

3.1.2  Integration  Of  The  Equations  Of  Motion 

The  methods  chosen  to  integrate  the  equations  of  motion  are  well  known, 
Reference  [3].  Two  methods  were  selected,  the  modified  Euler  and  the 
classical  fourth  order  Runge-Kutta  method.  They  were  chosen  because  of  then- 
stability  and  because  they  are  one-step  methods. 


Reference  [ 3 J ,  Bellman,  Richard,  Introduction  to  Matrix  Analysis.  New  York: 
McGraw  Hill  Book  Company,  1960. 
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One-stop  methods  have  the  advantage  that  only  information  within  the  interval 
of  integration  is  used  and  the  function  and  its  derivatives,  up  to  the  order 
of  the  method,  need  to  be  continuous  within  the  interval.  Multi-step  methods 
require  information  from  several  steps  and  the  function  and  its  derivatives, 
up  to  the  order  of  the  method,  need  be  continuous  over  all  the  intervals  used. 
Unfortunately,  the  later  condition  is  not  met  in  powered  flight  because  the 
acce 1 erat i ons  are  discontinuous  from  one  minor  cycle  to  another. 

Ihe  equations  of  motion  were  formulated  in  two  different  ways  and  the  two 
methods  of  integration  were  used  on  both  f ormulations. 

The  modified  Euler  method  is  as  follows 

Given  the  differential  equations 

y  --  f(y,  t.)  (3.5) 

y<V  '  yo 

where  y  is  the  unknown  variable,  in  general  a  vector,  and  specified  at  some 
initial  time,  then  by  the  modified  Euler  method 

“n.l  yn  '  h'(V 

Vl  -  y,  •  pIXP,,.,.  t„.,>  •  <<y„.  V> 

where  n  ’s  the  step  size  of  the  integration  Pn+j  considered  to  be  the 
predicted  vaiue  of  the  function  y  by  t  ,  y  ,  is  the  corrected  value  iiin 

(3.6)  it  is  _ietr  that  the  values  of  y  at  one  step  are  completely  determined 

by  values  at  the  reeding  step.  This  method  is  said  to  be  a  second  order 
-etnod  because  t  were  a  polynominal  in  t  of  no  higher  degree  than  two,  thru 

the  modified  Ei  method  would  integrate  the  differential  equations  exact :> 

if  the  arithmet  *"e  done  in  infinite  precision. 

any  function  y  which  lias  a  third  derivative,  the  truncation  error  incurred 
uy  using  this  method  is  given  by 
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where  t;  lies  somewhere  within  the  interval  of  integration. 


The  second  method  used  was  the  classical  fourth  order  Runge-Kutta  method.  It 
is  as  follows: 

Given  the  equations  of  (3.5),  the  solution  is  determined  by 

yn+l  yn  +  l[kl  +  k2  +  k3  +  k4J  (38) 

where 

kl  =  hf(yn*  tn) 

k2  =  hf(yn  *  **ki«  ln  +  ^ 
k3  =  hf(yn  +  hk2,  tp  ♦  *sh) 

k4  =  *  k3 '  ln  *  h)  • 

This  method  is  fourth  order  because  if  y  were  a  polynomial  of  degree  no  higher 
than  four,  then  (3.8)  would  give  the  answer  exactly  at  any  time  if  the  arithmetic 
were  carried  out  in  infinite  precision. 

Its  truncation  error  is  difficult  to  determine,  which  is  one  of  the  weaknesses 
of  the  method;  however,  for  any  function  which  has  a  fifth  derivative,  there 
exists  a  positive  constant  c,  such  that 

f.T  <  ch5  (3.9) 

where  h  is  the  step  size. 

With  the  truncation  error,  or  the  upper  bound  for  it,  expressed  as  it  is,  it 
is  possible  to  formulate  an  approach  which  will  determine  an  upper  bound  for 
the  total  error  incurred  in  integrating  the  trajectory.  Figure  1  shows  that 


J 


if  the  step  size  is  taken  large  enough,  the  truncation  error  is  much  larger 
than  the  round-off  error.  In  this  case  t:  ~  The  approach  was  to  run  the 

program  at  such  a  step  size  to  a  fixed  time  and  then  cut  the  size  in  half  and 
rerun  the  program  to  the  same  time.  By  subtracting  the  two  results  a  very  good 
estimate  was  obtained  of  the  error  incurred  by  using  the  larger  step  size. 

The  equations  are  as  follows: 

*c(t,  h)  =  xy(t)  +  eT(t,  h)  +  rR(t,  h) 

The  computed  value  at  time  t  for  a  given  step  size  h  is  equal  to  the  true  value 
Xy(t)  plus  the  truncation  error  and  round-off  error  for  the  step  size.  However, 
we  have  assumed  that,  at  the  step  size  chosen,  the  round-off  error  is  much 
smaller  than  the  truncation  error.  *;R(t,  h)  <  <  ;^(t,  h).  Tnerefore, 

xc(t,  h)  ~  xy(t)  ♦  tT(t,  h)  .  (3.9) 

Now,  in  both  the  modified  Euler  and  in  the  Runge-Kutta  methods,  the  truncation 
error  is  expressed  as 

cT(t,  h).  =  c.hk 

where  k  =  3  for  the  modified  Euler  method,  and  5  for  Runge-Kutta  method.  The 
subscript  i  indicates  that  the  error  term  applies  only  to  making  a  single  step. 
However,  when  the  step  size  is  halved,  then 

,t(t,  h/2)  j  =  '-i  hk 

assuming  the  ktr  4e-ivatives  are  nearly  constant  in  the  interval.  Since  two 
steps  need  to  be  t  »en  to  integrate  to  the  same  place  as  in  the  large  step, 

;.ne  truncation  e>ror  for  the  smaller  step  at  time  t  is  approximated  by 


*  T(t .  h/2). 


h 


k 
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(3.10) 


“  2k-l  £T(t’  h)i  ' 

Now  the  computed  value  at  half  the  step  size  is  equal  to 

xc(t,  h/2)  ~  xT(t)  +  eT(t,  h/2)  .  (3.11) 

By  subtracting  (3.11)  from  (3.9)  it  follows 

xc(t,  h)  -  xc(t,  h/2)  ~  &T(t,  h)  -  r-T(t,  h/2)  . 

Since  from  (3.10)  it  was  shown  that  the  truncation  error  at  each  step  for  the 
half  step  integration  is  a  fraction  of  the  error  at  the  full  step,  the  total 
truncation  error  up  to  time  t  is  mainly  that  from  the  larger  step  size. 
Therefore 


xc(t,  h)  -  xc(t,  h/2)  ~  cT(t,  h)  .  (3.12) 

The  above  approach  is  only  valid  when  the  truncation  error  needs  to  be 

determined  to  an  order  of  magnitude.  It  is  not  possible  to  make  a  precise 
determination  in  this  manner.  It  is  also  clear  that  the  assumptions  are  better 
met  by  the  higher  order  Runge-Kutta  method. 

For  the  above  approach  to  give  meaningful  results,  the  size  of  the  step  and 
the  half-step  must  be  large  enough  so  that  the  truncation  error  is  the  dominant 
term.  This  fact  is  verified  if,  for  a  smaller  value  of  h  and  half  its  value, 
closer  agreement  is  obtain  in  the  computed  values.  If  the  round-off  is  becoming 
the  larger  term,  for  smaller  values  of  h  the  differences  should  be  about  the  same 
or  become  larger. 

It  is  clear  from  comparing  the  number  of  times  the  right  hand  side  of  the  dif¬ 
ferential  equation  has  to  be  evaluated  in  (3.6)  and  (3.8)  that  the  Runge-Kutta 

method  will  take  longer  than  the  modified  Euler  method.  The  amount  of  work  is 

doubled  for  the  Runge-Kutta;  however,  because  the  Runge-Kutta  permits  larger 
step  sizes,  this  method  turns  out  to  be  the  most  economical  in  obtaining  the 
high  accuracy  desired. 
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3.1.3  Spline  Polynomials  For  The  Trajectory 

The  error  caused  by  interpolation  needs  to  be  controlled  as  carefully  as  the 
error  from  numerical  integration.  The  manner  for  determining  the  error  is 
rather  straight  forward.  After  the  coefficients  of  the  interpolating  polynomial 
have  been  determined,  it  is  tested  by  integrating  over  the  interval  of  inter¬ 
polation  with  a  very  short  step  size.  The  integrated  values  are,  consequently, 
very  accurate.  The  interpolated  values  from  the  polynomial  are  compared  to 
them. 

The  algorithms  for  determining  the  coefficients  of  the  polynomials  and  the 
programs  to  do  this  are  given  in  Reference  [1].  For  the  trajectory,  it  is 
possible  to  use  both  cubic  and  quintic  splines.  The  cubic  spline  is  determined 
by  using  position  and  velocity  at  each  end  of  the  interpolating  interval.  The 
quintic  spline  is  determined  by  using  position,  velocity  and  acceleration  at 
each  end  of  the  interval. 

The  polynomials  are  normalized  to  interpolate  over  in  interval  (0,  1).  This 
improves  the  round-off  error  appreciably.  With  the  position,  velocity  and 
acceleration  at  each  end  of  the  interval  the  normalized  polynomial  coefficients 
are  determined  in  the  following  manner. 

Let 


Reference  [1],  Thompson,  G.  T.  ,  Computation  of  State  Vector  and  Transition 
Matrix  for  TRAM.  Vandenberg  Air  force  Base,  Federal  Electric  Corporation,  WTR 
Division,  1977. 
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wne^e  p  ,  v  and  a  are  the  position,  velocity  and  acceleration  at  the  beginning 
of  the  interval,  and  p^,  v^  and  a^  are  position,  velocity  and  acceleration  at 
the  end  of  the  interval.  Each  of  these  quantities  is  a  vector  with  three 
elements,  h  is  the  length  of  the  interval,  h  =  t^  -  tQ. 

Then  the  coefficients  for  the  cubic  polynomial  are  determined  as  follows: 


c3  =  vi  +  v0  '  2(pl  "  po) 


(3.13) 


c2  "  pl  •  Po  ‘  vo  •  c3 


C1  =  vo 


co  =  po 


The  polynomial  is  then 


p(t)  =  cQ  +  Cjt  +  c2t2  +  c3t3 


and  is  computed  for  each  of  the  three  components  of  the  vector. 


The  quintic  polynomial  coefficients  are  determined  by  the  following  algorithm. 


a,  =  Pi  '  P  '  v  '  a/2 
1  1  ro  o  o 


(3. 14) 


a2  =  V1  '  vo  ‘  ao 


a3  =  (al  ’  ao)/2 


d^  =  -  3a^  +  6a^ 


ib 


^2  “  ~ 


al  '  d4  d5 


V2 


d,  =  v 
1  o 


d  =  p 
o  Ko 


The  polynomial  is 


5<t>  =  do  .  djt  .  d/  .  d3t3  .  d4t4  *  d6t5 


Since  the  polynomials  are  normalized,  position  and  velocity  at  any  time 
t,  tQ  <  t  <  tj,  are  obtained  from  the  cubic  by  these  equations 


p(t)  =  p( u( t ) ) 


(3.15) 


where 


v(t)  -  p'(u(t))/h 


0  u(t)  ~  (t.  -  t  )/h  <  1  . 


When  the  quintic  spline  is  used,  the  position  and  velocity  for  any  time  t  is 
given  by 


p(t)  -  ,(u(t)) 


(3.16) 


where 


v(t)  q 1 ( u( t) )/h 


0  <  u(t)  =  (t.  -  t  )/h  i  1  . 
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3.2.1  Mathematical  Formulation  Of  The  Transition  Matrices 

The  inatnematical  treatment  of  the  transition  matrices  is  given  in  Reference  [2], 
The  transition  matrices  are  related  to  the  fundamental  matrices  by  the  equation 

*(t.  s)  -  mif^s)  (3.17) 

and  the  fundamental  matrices  satisfy  the  differential  equation 

«i*( t)  =  F(t)f(t)  t'  <-  t  <  t"  (3.  18) 

t  * )  =  I 

where  I  is  the  identity  matrix  and  F(t)  is  given  by  (2.2). 

The  inverse  of  'f(t),  ^(t)  also  satisfies  a  differential  equation 

r^t)  =  -f'1(t)F(t) 

4i_1(0)  =  I  . 

t  ^(t)  may  be  found  by  integrating  the  differential  equations  or  by  finding 
Y(t)  and  inverting  it. 

In  this  study  only  those  terms  of  the  transition  matrix  which  are  associated  with 
position  and  velocity  are  considered.  This  means  the  F(t)  above  is  equal  to 
F  (t)  and  'F(t)  is  equal  to  4'  of  Section  6.3.  The  subscripts  are  dropped 

33  d,  <3 

to  simplify  notation. 

It  should  be  pointed  out  that  this  is  apparently  the  simplest  part  of  the  tran¬ 
sition  matrix  to  compute.  ^  of  Section  6.3  evidently  presents  some  significant 
problems.  Flow  these  problems  are  handled  will  appear  in  another  report. 


Reference  [2],  Brooks,  R.  A.,  Trajectory  Reconstruction  and  Analysis  Methodology, 
Vandenberg  Air  Force  Base,  Performance  Analysis  Department,  Federal  Electric 
Corporation,  WTR  Division,  1978. 
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3.2.2 


Integration  Of  I  he  Iranst.ion  Matrices 


The  problem  is  how  to  integrate  (3.18)  accurately  and  efficiently  and  then  to 
explore  what  accuracies  can  be  realized  with  the  chosen  methods.  It  is  important 
to  note  that  most  of  the  integration  time  will  be  spent  in  the  integration  of 
fundamental  matrices  simply  because  of  their  size.  For  whatever  methods  are 
chosen,  a  standard  with  which  to  compare  them  is  necessary  in  order  to  maintain 
the  desired  accuracy  I  he  inverse  of  l’(  t  )  will  be  integrated  the  same  way  T'(t) 
i  s . 

The  integration  of  the  fundamental  matrices  for  free  fall  was  examined, 

Reference  [1],  The  elements  of  the  matrices  are  t.he  same.  The  only  difference 
is  the  more  dynamic  and  discontinuous  powered  ‘light,  trajectory.  From  the 
mathematics,  it  appears  these  differences  would  make  only  second  order  changes. 

The  equations  to  show  this  follow 

t )  ^  F(tWt) 

t(0)  -  ; 

Where  I  is  a  identiy  matrix.  Whatever-  integrator  is  chosen,  how  well  it 

performs  depend*,  on  now  we1!  the  higher  derivatives  of  T  with  respect  to  t  oehave. 

’**<  t)  -  f  ( t  )■}■{ ;  i  *  f  ( t.  )i  ( t ; 

-  F(tMt)  ♦  f  (t)F(t)t(t)  . 

The  F ( t )*p ( t )  term  gives  a  clue  as  to  what  might  be  expected  from  the  integrators. 

F(t)  is  a  mat-  •  *  e  ijth  term  can  be  denoted  as  q.j(t);  therefore 

F(t)  ( <1  j  j  (  t )  )  i  .  J  i  ,  .  .  6  . 

From  the  definit  n  of  F(t)  it  is  clear  that  it  is  an  explicit  function  of  the 
position  and  velocity,  which  are  functions  of  time.  As  a  matter  of  fact  t  does 
no-  appear  explicitly  in  F(t);  therefore,  it  may  be  written  as 

Reference  [1],  Thompson,  G.  T.,  Computation  of  State  Vector  ana  Transition 
Matrix  for  TRAM,  Vandenberg  Air  Force  Base,  Federal  Electric  Corporation, 

WTR  Division,  1977. 
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F(t)  =  (gi j(x(t)))  i ,  j  -  1 ,  •  -  -  ,  6 


Then 


F(t)  =  (gij(x(t))) 


But 


The  term  x(t)  was  discussed  in  Section  3.1.2,  where  it  was  noted  that  the  powered 
flight  trajectory  had  greater  dynamics  and  that  the  Ap(t)  term  in  these  dif¬ 
ferential  equations  was  discontinuous.  In  this  expression,  however,  each  term 
in  x(t)  is  multipled  by  coefficients  that  are  either  zero  or  very  small.  For 

i  =  1,  2,  3  the  coefficients  are  identically  zero,  and  for  i  =  4,  5,  6  the  terms 
-12 

are  bounded  by  10  in  absolute  value.  As  a  result  the  greater  dynamics  and 
the  discontinuities  in  x(t)  will  have  a  second  order  effect  on  the  integration 
of  the  characteristic  equation. 

For  the  term  F(t)4'(t),  with  the  use  of  the  Schwartz  inequality 


llF(tWt)  ||  <  || 


11  x(t)  ||  ||  noil 


where  ||  •  j|  indicates  the  norm,  Reference  [3],  but  since 


II 


< 


||  x(t)  ||  <  10 

lint)  I!  -  i 


10 


-12 


5 


Reference  [3],  Bellman,  Richard,  Introduction  to  Matrix  Analysis,  pages  161  and 
162  ,  New  York:  McGraw  Hill  Book  Company,  1960. 
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then 

||  F(tWt)  !|  •  10'12  •  105  =  IQ'7  . 

Similarly,  we  can  see  by  inspection  that  higher  partials  of  g..(t)  will  be 
smaller  still.  The  higher  derivatives  of  n>(t)  are  discontinuous  but  the  dis¬ 
continuities  are  multiplied  by  very  small  partials,  and  as  a  result,  it  is 
expected  that  rather  large  step  sizes  might  be  used  in  the  integration  of  t(t). 

The  second  term  in  T(t)  is  F ( t )T ( t JT( t ) . 

Since  ||  F  ||  '  10  J ,  ■ !  H*(  t }  1 1  _  ],  and  this  term  is  continuous,  it  presents  no 

difficulty  in  the  integration  of  H'(t). 

Tne  methods  of  integration  chosen  were  the  fourth  order  Adams-Moul ton  and  Runge- 
Kutta  methods.  The  Runge-Kutta  does  not  differ  concept.ional iy  from  the  way  in 
which  it  was  discussed  in  Section  3.1.2,  although  the  differential  equations 
are  different  and  the  application  is  somewhat  different  The  explicit  algorithm 
is  given  by  (3.8). 


The  Adams-Moul ton  method  is  discussed  in  detail  in  Reference  [4],  By  it  the 
solution  of  (3.b)  is  as  follows: 

Vi  =  yn*  2-4  (55f(v  V  -  M,(Vr  Vi> 

4  37f<v2’  tn-2)  -  9,(Vr  v3)} 

•'rM  "  Vn  +  h  (9f(ynrl'  W  4  i<Jf<V  tn) 

St  ( y  . ,  t  . )  ♦  f(y  t  ,)) 

•'n-i  n-1  7n-2  n-2 

where  h  is  the  egration  step  size  and  t(y  t(i>  is  the  value  of  the  derivative 
of  y  at  t  . 


xeference  [4],  Lapidas,  Leon  and  seinfiled,  John  H. ,  Numerical  Solution  of 
Ordinary  Differential  Equations.  New  York:  Academic  Press,  Inc.,  1971. 
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Several  past  values  of  the  function  are  used,  and  if  the  method  is  to  be  effective, 

it  must  nave  high  order  derivates  that  are  continuous  in  the  interval  from  tn_^ 

to  tn+^.  In  integrating  from  t^  to  tn+^  the  method  first  computes  a  predicted 

value  y  ..  at  t  ,  and  then  the  derivative,  evaluated  at  t  .,  with  ,  is 
■'n+1  n+1  n+1  ■'n+1 

used  to  compute  the  corrected  value  yn  +  j  at  tn+^. 

Since  to  step  from  t  to  tn+^  only  two  evaluations  of  the  differential  equations 
are  required,  this  method  is  nearly  twice  as  fast  as  the  Runge-Kutta  method 
which  requires  four. 


This  method  is  known  to  be  very  stable  and  the  truncation  error  for  the  predictor 
and  corrector  is  given  by 


T 

P 


251  5  [5] 
720  y 


T 

c 
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.5  [5] 
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The  corrector,  as  is  always  the  case,  is  more  accurate.  Both  error  terms  depend 
upon  the  fifth  power  of  the  step  size  and  the  fifth  derivative  of  the  function 
being  integrated. 


The  first  three  values  of  the  variable  after  the  initial  conditions  have  to  be 
obtained  by  other  means  before  this  method  can  be  used  for  the  first  time. 
These  initial  three  values  were  obtained  by  the  Runge-Kutta  method. 


3.2.3  Interpolators  For  The  Transition  Matrices 

Two  interpolators  were  used.  One  was  the  cubic  spline  described  in  Section 
3.1.3.  The  other  was  a  interpolator  which  uses  the  function  and  the  derivative 
at  three  points.  This  is  not  a  spline  interpolator  in  the  correct  sense  but 
it  is  similar  to  one.  The  quintic  spline  described  in  Section  3.1.3  was  not 
used  because  the  second  derivative  of  the  fundamental  matrices  are  needed  for 
it  and  these  are  difficult  and  very  costly  to  compute. 
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The  cubic  spline  has  already  been  discussed,  but  the  quintic  interpolator  needs 
to  be  expressed  explicitly.  What  is  given  are  three  values  of  a  function  and 
its  derivative  at  equally  spaced  points.  So  then  these  values  are  given 

tup,  t'Ctp,  t(tp,  t'(tp.  tUp ,  'f'Up,  h.  h  =  t3  -  tj  . 

Compute  first 

rx  =  -  -  t'ftp  -  t(tp  r  t(tp 

r?  -  -  ht'Up  ♦  ht(t2) 

r3  =  -  ht'Up  -  t(t2)  ♦  t(t3) 

r4  =  -  ht'(tl)  ♦  ht(tp 

Then  compute 

ao  =  t(tp 

a}  =  ht'Up 

=  16  p  -  8»2  +  7r3  -  r^ 

a3  =  *  r4  *  4a2 

a^  -  32p  -  r  ^  -  3a3  -  7a^ 

r  3  -  -  a^  - 

"hen  t  te  po  1  yrn"  i 

2  3  4  6 

p(u)  s  *  a.u  *  a„u  +  a.u  ♦  a„u  +  a,u  ,  0  <  u  ^  1 

u  i  c  J  e  b  - 

u  =  (t  -  tp/h 


2? 


s  .»  l.o  lynomial  ,  normalized  over  the  interval  t^  to  t^  which  .matches  the  giver, 
functional  values  ana  derivatives  at  each  end  of  the  interval  and  at  the  mid-point 
That  is, 


P(0) 


'l'(t1) 


p'(0)/h  -  r(t.) 


P(>»)  ~  l'(t2) 


p'  (H)/h  -  r(t2) 


Pd)  =  n  t3) 


p'(i)/n  -  r(t3)  , 


'o  fim  the  interpoiated  value  of  any  element  of  the  fundamental  matrix  or  it 
derivative  at  any  value  t,  t^  t  <  t^,  u  is  computed  first 


u  =  (t  -  tp/h 


anu  then  p(u)  and  p'(u)/h. 


*Kt)  =  p  ( u ) 
T'(t)  =  p'(u)/h 


4.0 


NUMERICAL  RESULTS 


4.1  Computation  Of  The  Trajectory 

4.1.1  Integration  Of  The  Trajectory 

With  the  equations  formulated  as  in  Set  I  and  as  in  Set  II  discussed  in  Section 
3.1.1,  both  sets  were  integrated  using  the  modified  Euler  and  the  Runge-Kutta 
methods  described  in  Section  3.1.2  to  200  seconds. 

The  only  parameter  varied  as  both  sets  of  equations  were  integrated  by  the  two 

methods  was  the  integration  step  size.  The  step  size  was  changed  in  such  a 

way  that  either  a  minimum  total  error  was  reached  (as  described  in  Section  2.0) 
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or  the  error  in  position  became  less  than  10  feet  and  simultaneously  the  error 
in  velocity  became  less  than  10  ^  feet/second.  In  the  numerical  experimenta¬ 
tion,  the  errors  fell  below  the  bounds  before  the  minimum  error  was  reached. 

Table  1  gives  the  results  from  the  various  integrations.  The  equations  are 
identified  as  Set  I  or  Set  II;  the  step  size  and  the  method  of  integration  are 
indicated.  The  error  incurred  is  given.  The  error  was  determined  in  the 
manner  described  in  Section  3.1.2.  Also  given  is  the  time  taken  to  do  the 
integration,  h  is  the  step  size.  e  is  the  error  in  position,  and  is  the 
error  in  velocity. 
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TABLE  1 

TRAJECTORY  INTEGRATION 


MODIFIED  EULER 

RUNGE-KUTTA 

n 

SET  1 

SET  2 

SET  1 

SET  2 

.015 

T  -  411 
c  =  1. 5xl0~4 

t  =  10  6 

T  =  419 
c  =  8xl0‘6 

P  -  7 

r.  =  9x10 

V 

.03 

T  =  218 

f.  =  6xl0'4 

P 

t  =  4x10  0 

V 

T  =  231 

c  =  5x10  4 

P  -t 

c  =  9x10 

V 

T  =  346 

i.  =  6xl0"6 

P 

i.  "  1x10 

V 

T  =  351 

t:  =  7xl0'6 

P  _  O 

I.  =  3x10 

V 

.06 

T  --  97 

&  ---  10_1 

p  -5 

r.  =  5x10  3 

V 

T  =  119 

£  -10’3 

<.P  =  2xl0'5 

V 

T  =  168 

f.  =  4x10 

P  -5 

i.  =  1.5x10 

V 

T  =  170 

i.  =  9xl0"4 

P 

t  -  5x10  D 

V 

.09 

T  =  94 

c  =  3x10* 3 

P  -5 

c  =  2x10 

V 

*  In  Table  1  T  has  units  of  seconds,  t:  has  units  of  feet,  i.  has  units  of 

p  v 

feet/second,  and  h  has  units  of  second. 

The  times  given  must  be  divided  by  three  to  give  the  time  of  the  integration 
at  the  specified  step  size.  This  is  because  the  step  size  was  divided  by  2  to 
integrate  a  value  that  was  used  as  a  standard  for  comparison.  The  integration 
time  at  half  the  step  size  plus  the  integration  time  at  the  full  step  size  is 
about  three  times  the  value  at  the  specified  step  size. 

The  cases  where  no  data  is  presented  in  the  table  were  not  run  because  either 
the  error  had  become  too  large  or  too  small. 

4.1.2  Trajectory  Interpolation 

The  algorithms  for  the  cubic  and  quintic  splines  used  to  interpolate  the  traj¬ 
ectory  a^e  given  in  Section  3.1.1.  The  equations  to  compute  the  coefficients 
of  the  polynomials  are  given  by  (3.13)  and  (3.14). 
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Both  the  cubic  and  quintic  splines  were  used,  but  the  results  were  so  similar 
that  it  was  not  meaningful  to  treat  them  separately.  There  was,  however,  a 
difference  in  the  results  when  the  equations  of  motion  were  formulated  dif¬ 
ferently,  and,  of  course,  the  results  changed  markedly  when  the  step  size  was 
changed. 

Table  ?.  gives  the  results  of  the  numerical  experimentation  on  spline  i  nterpo  i  at  .  on 
It  gives  the  error  in  interpolation  obtained  by  subtracting  the  very  accurately 
integrated  values  from  spline  interpolation  The  parameters  are  the  same  as  in 
Table  1. 

TABLE  2 

SPLINE  INTERPOLATION 


SET  1 


T  =  946 

t  =  10'8 

P  -7 

t.  -10 
v 

T  =  506 


t:  =  10 

v 


SET  2 


».  =  10 

v 

T  =  519 
<  -  10' 


t.  -  10 

v 


T  =  249 


*  =10 
v 

T  =  136 

c  =  3xl0‘? 

P 

(  -  2 . 5x10 


T  =  73 

(.  =  3x10 

P  -5 

t.  =  10 
v 

T  =  39 

c  =  10'4 

P  -5 

c  =10 
v 


•«  2  Computation  Of  1  he  Characteristic  Matrices 

421  Integration  Ot  The  Characteristic  Matrices 

The  characteristic  matrices  and  their  inverses  were  integrated  in  two  ways 
One  was  by  the  Runge-Kutta  method  and  the  other  was  by  the  Adams-Mou I  ton  method 
which  were  described  in  Sections  i  1  2  and  3  2.2  Only  tne  character > st  ic 
matrices  associated  with  the  dynamic  terms  (t  )  were  integrated.  Reference  '2 j 

dd 

The  inverses  were  calculated  both  by  integration  and  by  inversion 

The  Adam»-Mou 1  ton  integration  was  allowed  to  run  over  the  er.t > r e  span  of  inte¬ 
gration  (100  seconds)  without  being  reinitialized,  whereas,  the  Runqe- Kuttu 
integration  was  reinitialized  to  the  identity  at  every  1  I ep 

The  standard  with  which  the  integration  ot  these  matrices  and  t he i •  'nverses 
by  both  methods  were  compared  to  was  the  result  of  Runge-Kutta  ’ntegrations 
witn  much  smaller  step  sizes. 

Since  the  point  of  view  is  different  in  the  two  different  methods,  the  tesa  ts 
of  the  integrations  are  presented  in  two  different  tables 

4.2  2  Interpolation  Of  The  Character) st ic  Matrices 

The  i nterpo i at i on  of  the  characteristic  matrices  was  accomplished  by  cubic  and 
qointic  polynomials  described  in  Section  3.2.2  When  the  Adams-Mou 1  ton  method 
was  useo  as  an  integrator,  both  the  cubic  spline  polynomial  and  qumtic  polynom 
was  used  The  standard  used  for  comparison  was  the  result  of  Runqe-Kutta 
inteyrat’on  with  small  step  sizes  The  polynomials  were  evaluated  at  the  step 
points  of  the  Runge-Kutta  integration 

When  the  Runge-Kutta  method  was  used  as  an  integrator,  only  the  cubic  spline 
polynomial  was  used  as  an  interpolator.  The  standard  for  comparison  was  the 
output  of  the  Runge-Kutta  method  small  with  step  sizes. 


Reference  [2],  Brooks,  R.  A  ,  T rajectory  Reconstruction  and  Analysis  Methodology 
Vandenberg  Air  force  Base,  Performance  Analysis  Department,  Federal  Electric 
Corporation,  WTR  Division,  1978. 
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4.2.  3 


Error  In  The  Characteristic  Matrx  Computation 


Tables  3  and  4  give  the  results  of  the  computation  of  characteristic  matrices 
for  the  Adams-Mou I  ton  and  Runge-Kutta  methods  respectively.  The  error  of  inter¬ 
polation  is  given  for  each  of  the  methods.  Both  the  cubic  and  guintic  poly¬ 
nomials  were  use  as  interpolators.  Since  the  results  were  nearly  the  same, 
only  the  values  from  the  cubic  were  considered. 

The  standard  matrices  were  subtracted  from  the  matrices  determined  by  inte¬ 
gration  and  by  interpolation  to  obtain  error  matrices.  The  norms  of  the  error 
matrice  were  taken  to  be  the  maximum  absolute  value  of  any  element  of  the 
matrix.  The  norms  of  the  error  matrices  are  given  in  the  tables. 

In  the  same  manner,  Reference  [1],  the  error  matrices  are  broken  up  into  fuu» 
parts  The  upper  left  nine  elements  are  those  errors  which  propagate  error-, 
in  pos’tion  into  position.  The  upper  right  nine  elements  submatrix  is  the 
matrix  which  propagates  errors  in  velocity  into  position  The  lower  left  nine 
element  submatrix  propagates  errors  in  position  into  velocity,  and  the  lower 
right  nine  element  submatrix  propagates  errors  in  velocity  into  velocity  The 
thirty-six  element  character i st i c  error  matrix  is  represented  by  four  number  . 
a  norm  for  each  of  the  submatrices  described  above. 


Re*e-»'  e  ,1],  Inompson.  G  T  ,  Computation  <>f  State  Vector  and  f rans  1 1 i on 
V  tor_TRAM,  Vandenberg  Air  Force  Base,  federal  Electric  CorporatTon,  toTR 
.i  i  .  TsTon  ,  19  7  7 
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TABU  1 

AOAMS-MOUt  JON  INI t  CHAT  JON 


ERROR  MATRICES 


STEP  i>IZt 

99 

1 . 98 

j.  9u 

INTEGRATION 

io"9  IO*8 

io'11  10-9 

io"8  io'7 

IO'10  IO'8 

IO'8  io'b 

iu'iG  lu'8 

INTERPOLATION 

io*9  io"8 

io"11  iu"9 

io'9  lu'  1 

10" lu  io'8 

_ 

io'8  io'b 

10* 10  10'8 

TAB. I  4 

RUNGE -KlM  ! A  iNiluKAl  ION 


E  RROR  MA I R I C  t  3 


3:1?  SIZE 

y9 

1 . 98 

.5. 9o 

IO'14  lu'19 

c- 

1 

O 

-  1  ?  -  0 

io  in  ' 

INTEGRATION 

IO'  18  10' 18 

10' 18  1G  14 

io' lii"1. 

INTERPOLATION 

IO' ^  ID' 11 

10- 19  10  M 

IO' 18  in-'1 

IO'14  lu' ‘ 1 

10-1“  U;-8 

I"  !  ''  lu"  11 

5.0 


CONCLUSIONS 


5.1  Conclusions  About  the  Trajectory  Computation 

The  trajectory  was  integrated  over  two  hundred  seconds.  From  Table  1  it  is 
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seen  that  the  only  times  that  an  error  in  position  of  less  than  10  feet  and 
simultaneously  an  error  in  velocity  of  less  than  10  feet/second  was  reached 
were  when  Equation  Set  1  and  Set  2  were  integrated  at  a  step  size  of  0.03 
seconds  with  the  Runge-Kutta  method  and  when  both  of  these  equation  sets  were 
integrated  at  a  step  size  of  0.015  seconds  with  modified  Euler  method.  These 
were  the  only  satisfactory  cases  of  all  the  cases  run. 

The  cases  chat  were  satisfactory  in  the  spline  interpolation  were  when  Set  1 
was  integrated  and  interpolated  over  an  interval  of  0.03  seconds  and  when  Set 
2  was  integrated  and  then  interpolated  over  intervals  of  0.03,  0.06  and  0.12 
seconds.  When  Set  2  was  interpolated  over  an  interval  of  0.24  seconds,  the 
position  accuracy  fell  well  within  the  requirement  but  the  velocity  error  was 
slightly  larger;  therefore,  this  interval  might  be  considered  satisfactory. 

5.2  Conclusions  About  the  Characteristic  Matrix  Computation 


Errors  in  the  computation  of  the  transition  vector  relate  errors  in  the  state 
vector  from  one  step  to  the  next.  If  we  bound  the  errors  in  position  and 
velocity,  then  the  total  error  can  be  bounded  by  multiplying  the  error  bound 
for  one  step  by  the  total  number  of  steps  taken. 


E 


T 


E  .N 
i 


step 


The  total  erroT  is  less  than  or  equal  to  the  bound  for  any  step  times  the 
total  number  of  steps  The  bound  for  step  i  is  given  by 


30 


where  E^  is  the  bound  for  the  computation  error  in  the  transition  matrix  for 
any  step  and  p£  and  v  are  the  bounds  for  the  errors  in  position  and  velocity 
for  any  step. 

With  the  same  approach  to  bounding  the  error  for  the  characters i t ic  matrix 
computation,  Reference  [1],  the  following  may  be  stated.  The  errors  in  position 
at  each  step  must  be  less  than  100  feet  and  in  velocity  less  than  one  foot  per 
second  and  since  there  are  less  than  1500  steps,  errors  in  the  characteristic 
matrix  computation  for  the  two  methods  of  integration  can  be  bounded. 

For  the  Adams-Moul ton  method  the  results  are  presented  in  Table  5. 


TABLE  5 

ADAMS-M0ULT0N  INTEGRATION  ERROR  FOR  TRANSITION  MATRICES 


tp  and  cv  is  the  bound  for  position  and  velocity  error  due  to  numerical  errors 
in  computing  the  transition  matrix. 

Since  the  interpolation  does  not  contribute  any  further  appreciable  error,  the 
values  in  Table  5  can  be  considered  the  total  error  for  this  method. 

From  the  above  tabled  results  the  following  conclusions  can  be  made.  The  Runge- 
Kutta  method  at  5.04  second  step  size  is  an  adequate  method  for  computing  the 
transition  matrix.  With  this  method  of  integration  the  cubic  spline  can  be 
used  sat i s  factor ly  to  interpolate  at  this  step  size. 

Reference  [1],  Thompson,  G.  T. ,  Computation  of  State  Vector  and  Transition 
Matrix  for  TRAM,  Vandenberg  Air  Force  Base,  Federal  Electric  Corporation,  WTR 
Division,  1977. 


This  step  size  is  in  agreement  with  the  analysis  which  anticipated  that  trie 
computation  of  the  transition  matrix  for  powered  flight  could  be  accomplished 
by  a  procedure  similar  to  that  used  in  free  fall. 


From  the  tables  it  is  apparent  that  the  interpolation  is  the  factor  which  limits 
the  step  size  for  the  Runge-Kutta  method 

As  for  the  Adams-Moulton  method,  it  is  clear  that  the  accuracies  obtained  by 
this  method  are  not  sufficiently  good  to  warrant  its  use.  It  must  be  remembered 
that  the  integration  was  made  over  100  seconds;  whereas,  in  the  Runge-Kutta 
the  integration  was  reinitialized  every  b  seconds.  Nevertheless,  even  when 
the  step  size  was  taken  to  be  99  seconds  the  accuracy  was  marginal.  This 
method  should  not  be  abandoned,  however,  because  it  has  the  potential  for 
integrating  as  well  as  the  Runge-Kutta  and  is  twice  as  fast. 

The  disadvantage  of  using  the  Adams  method  is  that  it  is  necessary  to  use  a 
different  starting  algorithm  for  the  first  three  steps  that  either  is  the  Runge- 
Kutta  method  or  is  as  costly  as  the  Runge-Kutta  method.  This  means  that  the 
integration  must  run  at  least  20  seconds  before  any  time  savings  are  made  at 
all,  when  a  five  second  step  size  is  used. 

For  the  Runge-Kutta  method  of  integration  we  have  a  similar  table 


TABLE  6 

RUNGE-KUTTA  INTEGRATION  ERROR  FOR  TRANSITION  MATRICES 


STEP  SIZE 

.99 

1.98 

3.96 

5.04 

c 

p 

-9 

°xl0 

2xl0'8 

2xl0~5 

2xl0*5 

i. 

v 

.  r 

. . 5x10 

,  1frio 

3x10 

2xl0'8 

3xl0"8 
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RECOMMENDATIONS 


Since  the  Runge-Kutta  error  is  less  and  is  faster  than  the  modified  Euler,  the 
recommended  algorithms  for  the  trajectory  are  the  Runge-Kutta  method  of  inte¬ 
gration  on  equation  Set  2  with  a  step  size  of  0.03  seconds  and  then  inter¬ 
polation  with  a  cubic  spline  for  no  more  than  0.24  second  intervals.  Because 
the  Runge-Kutta  error  is  less  and  is  faster,  the  recommended  algorithms  to 
compute  the  characteristic  matrices  are  the  Runge-Kutta  method  for  integration 
with  a  step  size  of  five  seconds  and  then  the  cubic  spline  interpolation  over 
the  same  interval. 
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